function [M, tList] = solveSteadyState(obj, tRange, npts)
    %solveSteadyState calculate M vector with periodic condition
    %   input:
    %       tRange: Specify a time range, which should contain complete periods of Omega. [in unit of s]
    %       npts:   Initial number of mesh points for bvp5c() solving. (Default value is 20.)
    %               If this value is too small, bvp5c() will increase it adaptively. Giving an "npts" that is
    %               either too small or too large will significantly slower the solving speed. The recommended 
    %               value is 10~30 points for each oscillation period.
    %   return:
    %       M: [3 x length(tList)] array, Mx = M(1, :); My = M(2, :); Mz = M(3, :);            
    %       tList: double row vector, time array [s]

    if nargin < 3
        npts = 20;
    end

    % check if tRange contains complete periods of all parameters
    tol = 1e-10;
    parList = {obj.parameters.T1,      'T1';
               obj.parameters.T2,      'T2';
               obj.parameters.OmegaX,  'OmegaX';
               obj.parameters.OmegaX,  'OmegaY';
               obj.parameters.OmegaX,  'OmegaZ';
               obj.parameters.Rp,      'Rp';
               obj.parameters.M0,      'M0';};
    for ii = 1:size(parList,1)
        err = obj.periodCheck(parList{ii,1}, tRange(1), tRange(end));
        if err > tol
            errMsg = sprintf('Residual returned by periodCheck() for "%s" is %.3g while the tolerance is %.3g.\n',parList{ii,2},err,tol);
            errMsg = sprintf('%sIt seems that the input "tRange" does not contain complete period of "%s".\n',errMsg,parList{ii,2});
            errMsg = sprintf('%sIt''s necessary that "tRange" contains whole periods of all parameters. Otherwise, the program will return a fake steady solution.',errMsg);
            error(errMsg); %#ok<SPERR>
        end
    end


    guess = @(t) [0.1*sin(2*pi*t); 0.1*cos(2*pi*t); 0.3];
    bcfcn = @(ma, mb) [ma(1)-mb(1);
                       ma(2)-mb(2);
                       ma(3)-mb(3)];
    BlochEquation = @(t, M) obj.blochEq(t, M, ...
                                        obj.parameters.T1, ...
                                        obj.parameters.T2, ...
                                        obj.parameters.OmegaX, ...
                                        obj.parameters.OmegaY, ...
                                        obj.parameters.OmegaZ, ...
                                        obj.parameters.Rp, ...
                                        obj.parameters.M0);



    tList0 = linspace(tRange(1), tRange(end), npts);                                
    solinit = bvpinit( tList0 , guess);
    options = bvpset('RelTol',1e-10, 'Stats','off', 'NMax', 2e4);
    bvpSol = bvp5c(BlochEquation, bcfcn, solinit, options);

    % refine tList for better plot appearance
%     tList = refineGrid_1D(bvpSol.x, 20);
    tList = linspace(tRange(1), tRange(end), length(bvpSol.x)*20);
    M = deval_fast(bvpSol,tList);


    obj.solution.type = 'SteadyState';
    obj.solution.tList = tList;
    obj.solution.mVector = M;
    obj.solution.sol = bvpSol;  
end